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Abstract 

Non-Markovian quantum state diffusion (NMQSD) is a non-relativistic but otherwise exact the- 
ory which expresses the reduced density matrix of an arbitrary subsystem, interacting hnearly with 
an uncoupled harmonic oscihator bath, as an average of diadics formed from state vectors which 
obey stochastic variational-differential equations. The vacuum radiation field can be represented 
as such an oscillator bath, and so this model is in widespread use in quantum optics. Prior to the 
development of NMQSD, exact subsystem solutions could only be obtained in a few special cases 
(e.g. spin-1/2, harmonic oscillator). Unfortunately, it has not yet been possible to obtain exact 
solutions to new problems using NMQSD due to the difficulty of solving the variational-differential 
equations. Here we show that these equations can be transformed into a pair of coupled nonlin- 
ear integrodifferential equations. We develop exact numerical methods for the integrodifferential 
equations and show that solutions can be readily obtained to good accuracy for quite general 
subsystems. We exactly solve various examples including tunneling in a double well representing 
molecular isomerization or racemization, suppression of fiuorescence from a two-level atom in a 
band gap, and intermittent fluorescence from a driven three level system representing electronic 
states of singly ionized magnesium. 



I. INTRODUCTION 



Few-level or few-body quantum systems such as single (or few chemically reacting) 
trapped atomic Q or molecular ions in ion traps, interacting with lasers and the vacuum, 
are systems of considerable current interest. Intermittent single ion fluorescence has been 

els of quantum computers^. These systems can all be accurately modeled as a subsystem 
interacting linearly with a bath of harmonic oscillators. Exact solutions for this model have 
only been obtained for restrictive cases where the subsystem consists of a spin-l/20,y] or a 
harmonic oscillator 0, [l^, and in a few other special cases. Theoretical studies of such sys- 
te.. have t.e.to.e a.:,,. .Ued o„ app.^^a. .a.te. e,ua.io,4 B , jump theoriesQ 
and Markovian stochastic Schrodinger equationsQ]. Nevertheless, a theory yielding exact 
solutions for this model for general subsystem Hamiltonians would have widespread applica- 
tions. Such a theory could also prove useful as a limiting case for development of theories for 
subsystem evolution under the influence of more general baths (e.g. single quantum-dot [l^ 
and single-molecule Q] fluorescence in condensed phase environments) for which theory is 
still in the early stages of development jl3 . 14 1, and for which intra-bath coupling is expected 
to play an important role [isj] . 

In a formal sense the problem of an arbitrary subsystem interacting linearly with an un- 
coupled oscillator bath has recently been solved. The first important contribution toward 
this theory showed that the influence functional has a particularly simple form. Subse- 
quently it was shown that the exact reduced subsystem density, expressed as a path integral 
weighted by the influence functional, can be stochastically unraveled as an average over 
diadics constructed from state vectors which obey a non-Markovian variational-differential 



wave equation 



17 



ISl |l9| . The resulting exact theory has been called non-Markovian quan- 



tum state diffusion (NMQSD), and from NMQSD one can in principle obtain exact solutions 
of the subsystem-oscillator bath model for any subsystem. 

In practice, no new exact solutions have been obtained using NMQSD. While exact 
evolution equations can be formulated for any subsystem, these equations contain variational 
derivatives with respect to the complex colored noises which emerge during the stochastic 
unraveling. Variational-differential equations (VDEs) have not received the same attention 
as partial differential equations and much less is known about their properties and solutions. 
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As a consequence the evolution equation has proved impossible to solve - even numerically 



except for the few models for which exact solutions were previously know n ll 7 . 



3. 



This 



impasse has led to the introduction of various perturbative expansions |2Q|,I 

I21I and numerical 

approximation schemes[22]. However, direct and exact methods for solution of the NMQSD 
equations remain an important goal. 

Direct solution of the NMQSD equations may eventually be possible, but for the present 
methods which circumvent the problem by eliminating the variational derivative seem most 
promising. This manuscript introduces a pair of coupled nonlinear stochastic integrodiffer- 
ential equations which we show are exactly equivalent to the stochastic VDE of NMQSD. 
Integrodifferential evolution equations can be converted to ordinary or partial differential 



equations 
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24| . Hence, the pair of stochastic integrodifferential equations can be converted 



to stochastic ordinary or partial differential equations which are easier to solve. Recently 
developed numerical methods allow efficient solution of stochastic (ordinary and partial) 
differential equations (SDEs)^, 2^ even for quite large systems of equations 22 1 . On the 
basis of the transformed equations and using these SDE integrators we develop exact numer- 
ical methods for solving NMQSD problems. Example calculations for which exact solutions 
are known are used to verify the accuracy and efficiency of the new methods. We also 
numerically solve three new problems which were previously intractable. 

Equivalent linear and nonlinear formulations of NMQSD exist. The nonlinear version 
of the theory preserves the norm of the state vector, which results in improved Monte 
Carlo convergence. Accordingly, we also consider linear and nonlinear reformulations of 
NMQSD. The simpler linear version is considered in section II, while the nonlinear case is 
treated in section HI. Simple example calculations are used to illustrate the two approaches. 
More interesting applications to a tunneling problem representing molecular isomerization 
or recemization, to suppression of fluorescence from a two-level atom in a dielectric band 
gap, and to a driven three level system exhibiting intermittent fluorescence, are considered 
in section IV. 
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II. LINEAR NMQSD EQUATIONS 



Consider for simplicity a subsystem-bath model with a single subsystem coupling operator 
L. The total Hamiltonian is 

Htot = H + gULal + L^a^) + ^ uala^ (1) 

where we are using units in which h = 1. Here H is the subsystem Hamiltonian, g^^ is a 
coupling constant for an oscillator mode of frequency u, and the rest of the notation should 
be clear. The generalization of NMQSD and our results to multiple subsystem coupling 
operators is straightforward, so we will confine our attention to this simplest case. 
In NMQSD the evolution of the state vector ipt is governed by the linear VDE 

^ = -tH + z,L - Lt /* ds ait, s) ^ (2) 
at Jo ozg 

where zt is a complex colored noise with correlation function 



a{t,s) = M[zlzs\. (3) 
In the limit of zero temperature s) = J2w 9'i^~^'^^^~'^^ (see Refs. 

01 



for the non- 



zero temperature formula). Here M[. . .] denotes the average over different realizations of 
the noise. The exact reduced density matrix pt of the subsystem is given as an average of 
diadics via pt = M[\tpt){tpt\]. ^ 

The solution ipt of Q is known to be an analytic functional of the noise Zt^lSj and so 
the variational derivative is well defined. The presence of this variational derivative does 
however make the evolution equation difficult to solve. In some simple cases one can guess 
the form of ^ and use a self-consistency requirement to find solutions of Eq. Q. Other 
than these few examples, it is generally unclear whether it is possible to directly solve the 
VDE (j2)) since numerical algorithms for VDEs have not yet been developed. Fortunately, 
it is possible to reformulate the theory in terms of more easily solved integrodifferential 
equations. 

To start with our reformulation, we introduce a non-unitary propagator U{t, s) which 
evolves a state vector from time s to time t. If we consider subsystem evolution from a state 
ipo at time then ipt = U(t, 0)ipo- From Eq. (j2)) we deduce that 

dUit,0) ^ , ^ ft^^ ^ 6Uit,0) 



dt 



iH U{t,0) + ZtLU{t,0)-L^ f ds a{t,s) , ^ (4) 

Jo OZs 



Next we need to eliminate the variational derivative ^^j^- In particular, we will show that 

6U{t,0) 



U{t,s) L U{s,0). (5) 



Consider for simplicity the case where L has a complete eigenbasis \x) which can be used 



to construct a path integral for the propagator along the lines followed in Ref. [17[. The 
path integral representation of U{t, 0) is a sum over paths weighted by an exponential whose 
argument includes a stochastic term Jq du ZuXu, where L\xu) = Xu\xu) and for each 
value of u denotes an element of some complete basis (i.e. the path integral is over all Xu 
for each value of u between and t)\v\. 18, Q|. The variational derivative of U{t,0) with 



respect to Zg thus brings down a prefactor Jq du inside the path integral. Using the 

fact that 

^ = 5(n - .) (6) 

dzs 

and removing the closure relation for Xg at time s we then obtain Eq. (0). 
Now Eq. and the semigroup property 

U{t,s) = U{t,0)U{s,oy^ (7) 

allow us to rewrite Eq. in the form 
dU{t,0) 



dt 



H U{t, 0) + ZtL U(t, 0) - U{t, 0) / ds a{t, s)U{s, Oy^L U{s, 0), (8) 

JO 



which is a nonlinear integrodifferential equation for U{t,0). This equation is exact and 
entirely equivalent to Eq. Q. No loss of generality is incurred in working with Eq. (jSI), but 
in a finite basis set the number of equations arising from (jHl) is the square of the number 
of equations arising from Eq. Moreover, Eq. (jSJ is clearly nonlinear whereas Eq. 

(j21) is linear in at least some cases. These apparent defects are unfortunate, but one must 
recognize that Eq. (jSJ is solvable while Q is not, and the computational costs of solving 
(jHl) are quite reasonable for a large class of potentially interesting problems. We solve three 
such examples in section IV. Thus, implementations and applications of (jH]) and its norm- 
preserving generalizations are the focus of this manuscript. 

Direct inversion of U{s,0), required for (jSJ, can be avoided at the expense of adding 
a second equation. Consider the dynamics of U{0,t) = U{t,0)^^. Using the fact that 
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U{0,t)U{t,0) = 1 and differentiating with respect to t gives 

dU{0, t) ^^^^ ^ ^^^^ t){~iHU{t, 0) + ZtLU{t, 0) 

(Jjv 

rt 



L^U{t,0) f ds a{t,s)U{s,0)-^LU{s,0)} = 0. (9) 
Jo 



Multiplying on the left by U{t, 0) ^ then gives 



dU{0,t) 



dt 
or 

dU{0,t) 



+ U{0, t){-iH + ztL - L^U{t, 0) /* ds a{t, s)U{s, 0)-^LU{s, 0)f/(0, t)} = (10) 

Jo 



tUiO,t)H - ZtUiO,t)L 
dt 

+ U{0,t)L^U{t,0) [ ds a{t,s)U{s,0)-^LU{s,0)U{(},t) (11) 
Jo 

which is again of integrodifferential form. Note that both equations (jHl) and (jlip involve 
only U{t, 0) and U{t, 0)"^ so that the pair is closed. 

Changing notation to Ut = U{t,0) and Uf-^ = U{0,t) we can rewrite equations ^ and 
(HD) as 



dUt 
dt 
dUf' 



-iH Ut + ZtL Ut - Ut t ds a{t, s)U^^L U, 

Jo 

iUt'^ H - ztUt'^L + Ut'^L^ Ut f ds a{t, s)U;^L Us Ut^ (12) 

Jo 



dt 

which is a closed set of integrodifferential equations. We will now show that Eqs. ()12|) can 
be transformed into sets of ordinary or partial differential equations. 

The most efficient set of transformed equations depends on the properties of the memory 
function. In section IIA we assume that the memory function consists of a few terms of 
exponential form, i.e., 

m 

a(t, s) = Y, ^ .g-7,|t-^|g-i^,(i-s) (13) 

where Aj and 7^ are positive numbers. The terms in Eq. ()13|) do not in general correspond to 
physical bath oscillator modes. Instead Eq. ( 2 ) can be viewed as a best fit to the memory 
function, obtained by nonlinear least squares 23 1 or other techniques In many cases the 



number of required terms m can be quite small. The general case where m is very large, or 
where a(t,s) cannot be written in the form (fTSj) . is considered in section IIB. To illustrate 
the application of these methods we numerically solve a number of example problems. For 



each example problem, Eqs. (1121) a re expressed in ordinary differential form and solved using 
stochastic integration methods j25l. l26|. 
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HI- 



We chose 



Colored noises can be generated using a variety of techniques |21l. 
memory functions for our examples which can be expressed via Eq. (jl3j) . The complex 
colored noise is then generated via Zt = J2JLi by integrating the stochastic differential 
equations 

' ■ (14) 



from — oo (in practice some large negative time) to time t = 0, and then from t = onward in 
combination with transformed versions of Eqs. p2|l . Here are complex Wiener processes 
satisfying M[dWrdW^] = Vdt6j,k5t,s and M[dWt^dW^] = 0. One can show that 



27,A, 



dWi e-^^(*-")e-*'^^(*-"\ 



(15) 



= J2Y=i^iy aiid Eq. yield the correct memory function ()13p in the mean. The 
stochastic differential equations (fT^ are obtained by differentiating expressions p5|) for 
Sets of stochastic differential equations like (jl4j) and the transformed versions of (|12p can 
be solved to any required tolerance using recently developed high order variable stepsize 



integration methods 



25 



A. Sum of exponentials 

In many cases it is possible to expand the memory function as a sum of a few exponentials 
of the form Eq. ()13|). Fits to expansions of this type can be obtained using nonlinear least 
squares algorithms. In practice, obtaining good fits with nonlinear least squares can be 
a frustratingly time consuming endeavor. The nonlinearity of the fitting function makes 
minimization algorithms highly sensitive to the parameter search domains. In such cases, it 
may be possible to obtain expansions of the same form using other techniquesj3]. 

Once this expansion is known we can define operators 

Vtj= f ds Aje-^^^'-'^e-'''^^'-'^U;'LUs (16) 

J 

such that the full set of equations in ordinary differential form becomes 

dlJ ™ 
dt ^tl 

7 



at .^^ 

^ = -(7, + ^^,)yt, + A,f/r^L/7,. (17) 

Efficient implementation of these equations requires that Y^=\ Vtj be computed ffist and 
stored temporarily. From this quantity UtJ2Y=i^t,j and L''UtJ2Y=i^t,j can be calculated. 
Numerically it is easier to use the fact that UfUf'^ = 1 and hence that dUf^/dt = 
—Uf^dUt/dtUf^ to find dU^^/dt than to program the above equation. 

To illustrate the utility of the method we will apply the transformed equations to a few 
problems where exact solutions are known. In all cases the SDEs were solved using the 



ANISE software 



2y|. We also used a single processor for each calculation, but of course one 



of the primary numerical advantages of quantum state diffusion theories is that they are 
inherently parallel and should ideally be implemented on clusters. With the exception of 
example IIC the calculations took only a few minutes. 



1. Example II A 

Consider the special case where H = {uo/2)az and L = Xa^ and pick s) = 
(7/2)e~'*'^*~*-' with uj = 1^ X'^ = 2uj and 7 = cj. We chose the initial condition \iI)q) = 
+ where <Jz\i) = {—iy~^\i) for i = 1,2 (this is the same convention as in Ref. 

ISf) . We calculated {l\pt\2) using Eqs. (fTTj) and an average over 10000 trajectories. In 
Fig. 1 we plot the real (solid curve) and imaginary (dotted curve) parts vs time against the 
corresponding known exact results[18] (dashed and dot-dashed, respectively). The accuracy 
is already very good. In Fig. 2 we plot the diagonal elements of pt which are supposed to be 
constant for this model. We see that the numerically calculated (l|pt|l) (solid curve) while 
initially equal to 5/7 (dashed line) decays away from this value as time proceeds. Similar 
results are observed for the element (2|p(|2) (dotted curve) which starts at 2/7 (dot-dashed 
line) and grows. Hence, for 10000 trajectories convergence is still incomplete for the diagonal 
elements. We will show in section III that introducing norm-preserving equations will lead 
to faster convergence. 

We also show in Fig. 3 the memory function calculated numerically using the stochastic 
process Zt (solid curve) and the exact memory function .5e~* (dashed curve). Agreement is 
satisfactory. 
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FIG. 3: a{t, 0) for Example IIA 




2. Example IIB 



Now consider the case where H = {uj/2)az and L = {X/2){ax — icTy) (Example IIIB in 
Ref. and again pick a{t,s) = (7/2)e~'^*^*^*^ with u = 1, X"^ = 2uj and 7 = cj. For 

the initial condition \iPq) = + |2)) we again calculated {l\pt\2) using Eqs. (fT7|) and 

an average over 10000 trajectories. In Fig. 4 we plot the real (solid curve) and imaginary 
(dotted curve) parts vs time against the corresponding known exact results |l8| (dashed and 
dot-dashed, respectively). Good agreement is obtained. Figure 5 shows the diagonal matrix 
elements (l|p(|l) and (2|p(|2) plotted against their corresponding exact results. The results 
here are also good. 



3. Example IIC 

Finally, consider a harmonic oscillator subsystem with H = ua'^a and L = Xa. We chose 
an initial condition \iPq) = ^^^\0) + ^|1) (where a^a\n) = n\n) for n = 0, 1, . . .) and picked 
a{t, s) = (7/2)6^^*^*^*^ with tu = 1, X = uj and 7 = u;. A basis set consisting of the first five 



10 




11 



oscillator states was employed. We computed 10000 trajectories to construct the average. 
The cpu time was about 114 minutes on a 600 MHz Alpha processor. In Fig. 6 we plot 
the real and imaginary parts of (l|pj|2) (solid and dotted) vs their corresponding exact real 
and imaginary parts (dashed and dot-dashed). In Fig. 7 we show that diagonal matrix 
elements (l|p(|l) and {2\pt\2) (sohd and dotted) plotted against their corresponding exact 
results (dashed and dot-dashed). The results are uniformly good. 




B. General memory function 



Here we consider the opposite case where the memory function cannot be efficiently 
represented as a sum of exponential terms like (jl3|) . In this case it may still be desirable to 



generate the colored noise using Eqs. 
spectral methods and other techniques 



21 



|_jbut otherwise the noise can be generated with 
The key point is that Eqs. (fT7j) will not 



be efficient and some new approach must be explored. 

Equations (fT^ are of integrodifferential form. Such equations can be solved exactly using 
a recently developed method j2J]. The trick is to convert the equations to partial differential 
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FIG. 7: {l\pt\l) and {2\pt\2) for Example IIC 




form. To do this we introduce a new variable u (where u G (—00, 00)) and define 

Vt,u = f ds a{t + u, s)U;^LUs. (18) 
Differentiating we find that the new set of equations has a partial differential form 

-iHUt + ZtLUt - L^UtVt,o 



dXh 
dt 

dt 
dVt,u 
dt 



a{t + u,t)U^^LUt + 



du 



(19) 



which is exactly equivalent to the original integrodifferential set p2|) . The partial deriva- 



tive can be evaluated using fast Fourier transform methods or using discrete variable 



du 



representations 



24 



0], or by employing a harmonic oscillator basis. In the first two cases 
the variable u is represented on a grid of points, while a discrete basis is used in the third, 
and so Eqs. (|19p revert to ordinary differential form and hence they can also be solved using 
the method of Ref. Ml or using the ANISE software Q- 

In practice there are advantages to modifying the above equations somewhat by instead 
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defining 

Vt,u = fiu) f ds a(t + u, s)U;'LUs (20) 
Jo 

where f(u)~ ^ where f{u) decays rapidly with u. These can lead to a smaller required 
basis setj2J| and hence faster algorithms. The modified Vt^u obeys 

= j{u)a{t + u,t)U^ LUt + 



-Vt,u 



(21) 



dt ' '"'"'^ ' ' du f{u, 

while the equations for Ut and Ut^ are unaltered. We will not explore the issue of which is the 



optimal form for /(m), but merely note that some success has been had with /(m) = e-""'|2J. 
In the case where the u degree of freedom is represented in an oscillator basis, f{u) should 
be chosen so that the oscillator matrix elements {n\f{u)a{t + u,t)\m) are nonzero only for 
the first few basis functions. 

Finally, we note that it is possible to use a similar trick to generate the complex noise Zt 



when it is stationary, i.e. when s) = c(t — s] 



21 



29|. Define 



f dWs R{t - s) 



(22) 



with R{t) = for t < 0, and where dWg is a differential Wiener process with properties 
M[dW;dWs] = Vdi6t,s and M[dWtdWs] = 0. It follows that 



from which one can show that 

Rit) 



dr R*{T)R{T-t) 



(23) 



/27r 



duj G{Lj)e 



iujt 



1 



(24) 
(25) 
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29|). 



One can then choose G{uj) = ^J\G{uj)\'^ (see Refs. 

Now that we have the representation (j^^ for the noise we can define 



zt,u = f{u) f dWs Rit + u-s) 

J —oo 



(26) 



with /(O) = 1 such that zt = Ztfl. Differentiating Eq. ()2(i|) we then obtain the partial 
differential equation 

^dzt^u f'{u) 



dztM 



du f{u) 



zt,u]dt + fiu)Riu)dWt 



(27) 
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which can be solved in concert with Eqs. ()19p and ()21|) . 

Obviously implementation of the rather sophisticated approach outlined in this section 
is quite a bit more involved than that for exponential type memory functions. However, all 
aspects of this treatment are exact and similar calculations have been shown effective for 



other types of integrodiffential equations 
approach to a subsequent manuscript. 



2J|. Hence, we will leave a detailed study of this 



III. NONLINEAR NMQSD EQUATIONS 

While we have obtained mostly good convergence and accuracy using the linear NMQSD 
equation, the theory can also be formulated in terms of a norm preserving nonlinear VDE. 
The introduction of norm-preserving equations is important since their solutions may have 
physical significance More practically, the norm-preserving equations in many instances 
yield faster convergence of the mean with the number of trajectories. The norms of the non- 
norm-preserving equations go to zero with probability one and so for long time dynamics, 
contributions to the mean tend to come from a few unusual trajectories potentially making 
convergence slow. The computational cost per trajectory for the norm-preserving equations 
is only a little greater than that of Eqs. (|T2|l . Moreover, the non- norm-preserving equations 
are themselves nonlinear and so there is no apparent drawback to the norm-preserving 
reformulation of NMQSD. 

The norm-preserving formulation of NMQSD is obtained via a two step process consisting 
of a Girsanov transformation followed by normalization of the wave vector. The details are 
given in Ref. j^. Defining Ut through ipt = Utipo, and substituting into the norm-preserving 
wave equation of Ref. Q], it can be deduced that 

^ = -tHUt + izt+ tdsa*it,s){L^)s) {L-{L)t)Ut 
at Jo 

- (Lt - (Lt)t) Ut f ds a{t,s)U:^LUs 
Jo 

+ OjoIUHl^ - {L^)t)Ut f ds a{t,s)U;'LUs\iJo) Ut (28) 

Jo 

where {L)t = {ipt\L\ipt)- Note that Ut now depends on the initial wave function ipo. As in 
the linear case there is also an equation for the inverse Uj~^, 



dUt-' 
dt 



iU;^H-{zt+ f ds a*it,s){L^)s) Ui' (L - {L)t) 
Jo 
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+ U-\L'^ - {L^)t) Ut f ds a{t,s)U;'LUs U^^ 

Jo 

- {tfJoluHL^ - {L^)t)Ut f ds a{t,s)U:'LUs\^o) U;\ (29) 

Jo 

These integrodifferential equations can again be re-expressed as sets of ordinary or partial 
differential equations. Again we consider two cases. 

A. Sum of exponentials 

If the memory function can be expressed as a sum of exponentials via p3|) then we can 
again define operators 

Vtj= f ds Aje-^^^'-'^e-'^^^*-'^U;^LUs (30) 

•J 

such that the full set of equations in ordinary differential form becomes 

dU ™ 

^ = -iHU, + {z, + Y.yt,3) (L - {L)t) u, 
at .^^ 

m 

i=i 

m 

+ {^l;o\uJ{L^-{L^)t)UtY.yt,J\^o) Ut 

i=i 

dU~^ ™ 

m 
m 

- {i^oWKl'^ - {L'^)t)UtY.yt,3\i^o) u^' 

^ = -{^,-iu:,)yt, + ML^)t (31) 

where yt,j = dsAje~^^^*~^^ e''^^^*~^^ {L^) s- A key to efficient numerical implementation of 
these equations is evaluation and temporary storage of Z^jLi from which Ut YJjLi 
and (Lt - (Lt)t) UtY.J=iVt,j, {^PolUJiL^ - {L^)t)UtET=iytJ^o) can be evaluated. Again 
dUt^/dt = —Ut^^dUt/dtUt^ should be used to find dU^^ /dt. The conservation law {■ipt\'ipt) = 
1 proves very useful for debugging code for Eqs. (jHH). When Tr{L} = an additional check 
can be made to verify that Tr{V^ j} = 0. 
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1. Example III A 



Consider again the case where H = {uj/2)az and L = \az and pick a{t, s) = ('j/2)e~"'^^~^^ 
with u = 1, = 2uj and 7 = cj. For the initial condition \iPq) = + ^|2) we 

calculated {l\pt\2) using Eqs. (|T7j) and an average over 10000 trajectories. In Fig. 8 we plot 
the real (solid curve) and imaginary (dotted curve) parts vs time against the corresponding 
known exact results 18] (dashed and dot-dashed, respectively). Once again we used the 

FIG. 8: (l|pt|2) for Example IIIA 
0.5 I 1 1 1 1 1 1 




-0.1 I 1 1 1 1 1 1 

0.5 1 1.5 2 2.5 3 

t 

ANISE software 26] to solve the equations. Norm was preserved to machine precision for 
individual trajectories. 

As in the case of the non-norm-preserving equations, convergence of the off-diagonal ma- 
trix element is good with 10000 trajectories. In Fig. 9 we show the diagonal matrix elements 
(l|p(|l) and (2 1 Pi 1 2). For the non- norm-preserving equations we obtained poor convergence 
for 10000 trajectories. However, for the norm-preserving equations the convergence is quite 
good. The non-preserving equations do indeed seem to lead to faster convergence as claimed 
in Ref. UM- 



17 



FIG. 9: {l\pt\l) and {2\pt\2) for Example IIIA 
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2. Example IIIB 

Now consider H = {uj/2)az and L = (\/2)(ax—icry) and again pick a{t, s) = (7/2)e~'^*-*~*-' 
with uj = 1, = 2u! and 7 = 0;. For the initial condition \ipo) = + |2)) we calculated 

(l|pf|2) using Eqs. (fT7|) and an average over 10000 trajectories. The real and imaginary 
parts (solid and dotted) are plotted in Fig. 10 against exact results (dashed and dot- 
dashed). Agreement is good. In Fig. 11 we show the numerically computed diagonal 
matrix elements (l|pt|l) and {2\pt\2) (solid and dotted) against exact results (dashed and 
dot-dashed). Convergence is again good. 

3. Example IIIC 

Now consider a harmonic oscillator subsystem with H = uoa^a and L = Xa. We chose an 
initial condition iV'o) = ^|1) + ^|2) and picked a{t, s) = /2)e^^^*-'^ with uj = l,\ = uj 
and 7 = u;. Convergence is rapid and so we used only 1000 trajectories. The cpu time was 
about 12 minutes on a 600 MHz Alpha processor. In Fig. 12 we plot the real and imaginary 
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parts (solid and dotted) of (l|pt|2) vs their corresponding exact real and imaginary parts 
(dashed and dot-dashed). In Fig. 13 we show the diagonal matrix elements (l|pt|l) and 
(2|p(|2) (solid and dotted) plotted against their corresponding exact results (dashed and 
dot-dashed). The results are good. 



FIG. 12: {l\pt\2) for Example IIIC 
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B. General memory function 

When the memory function cannot be represented as a sum of exponentials of the previous 
form, it is again possible to rewrite the equations in partial differential form. We again 
introduce a new variable u and define 

Vt,u = r ds a{t + u, s)U-^LUs (32) 
Jo 

such that the new set of equations has partial differential form 

^ = -tHUt + {zt + yt,o) {L - {L)t) Ut 
- (Lt - (Lt),) Ut Vt,o 
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FIG. 13: and (2|pt|2) for Example IIIC 




where 



dt 



+ {'UUliL^^ - {L'^)t)UtVtM Ut 



- a{t + u,t)Ui'LUt+^^'' 



t.u 



dt 
dyt,u 
dt 



du 

dyt,u 
du 



(33) 



yt,u= f ds a*{t + u,s){L''),. (34) 



The same considerations regarding adding a damping factor and generation of the noise Zt, 
discussed in section IIB, apply here without modification. 



IV. NEW EXAMPLES 

We now consider three new problems which it was previously impossible to explore exactly 
using NMQSD. Our treatment of these examples is not meant to be exhaustive. We simply 
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wish to show that interesting issues can be explored using the computational implementation 
of NMQSD introduced in sections II and III. We choose to use the norm-preserving version 
of the theory. This improves convergence and also allows us to study individual trajectories, 
which may have some physical significance Specifically, it has been shown that NMQSD 
can be interpreted as a hidden variable theory in which the noise Zt represents a hidden- 
variable of the bathjs^. We shall see that individual trajectories do indeed behave in 
ways consistent with physical intuition. Whether their non-quantum-mechanical statistical 
properties (e.g. M[|(-?/'o|V^i)P] = (V^oIaIV'o) can be predicted with standard quantum theory, 
while M [I ('?/'o|'^t) 1^] "~ ^[|(V'o|V'<)P]^ cannot) are in agreement with experiment remains an 
open and interesting question. 



A. Example IVA 

Here we consider the dynamics of a symmetric double well representing a reaction co- 
ordinate of an isomeric or chiral molecule, interacting with the radiation field. Of course 
realistic chemical environments contain sources of decoherence other than the radiation field, 
but the example is still of interest. When an ensemble of such systems is prepared in an 
initial achiral state, interaction with the radiation field is expected to drive the population 
to a symmetric final distribution. Indeed, clocks for dating amino acids have been proposed 
on this basisQ]. When the barrier height is low (e.g. in NHDT) individual molecules are 
observed in states which are superpositions of left and right handed states (e.g. the ground 
state of the double well). When the barrier height is large, individual molecules are found 
in left handed states or right handed states but not normally in superpositions (although 
superpositions can in principle be prepared js^). This is unusual because the eigenstates 
of the double well are superpositions of left and right handed states. Environment induced 
superselection rules are sometimes invoked to explain this effect j^, 0] . We will now explore 
the predictions of NMQSD in these two cases. 

Consider a quartic oscillator subsystem with H = uj[a^a — {3/8){a^ + a)^ + e{a^ + a)^] and 
L = Xa. This Hamiltonian corresponds to a symmetric double well potential and so could 
represent a reaction coordinate for isomerization or racemization of a molecule. We choose 
an initial condition = :^(|0) — |1)) which means the particle is initially localized in the 
left well, and pick a{t, s) = (7/2)6^^'^*^'*^ with u = 1, e = hcu/Eh, \ = cu and 7 = where Ef, 
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is the activation energy of the barrier. The parameter e controls the effective barrier height. 

First consider the case where e = .692, which is typical for proton transfer 
isomerizationjs^l. In Fig. 14 we plot the probability density {x\pt\x) against x (in units 



of ^jTi/muj) for times t = (solid), 2 (dashed), 4 (dotted), 6 (dot-dashed) and 12 (double- 
dashed). Dissipation and decoherence drive the population from a nearly pure left-handed 
state to a symmetric mixture. The calculation was performed in a basis set of the lowest 
five harmonic oscillator states. A total of 1000 trajectories were included in the average. 
The calculation again took about ten minutes. In Fig. 15 we plot {x\pt\x) against x, com- 



FIG. 14: Relaxation of {x\pt\x) vs x for Example IVA with e = .692 
0.7 




puted for 1000 and 10000 trajectories, at times t = 2 (solid and dashed), 4 (dotted and 
dot-dashed) and 12 (double-dashed and triple-dashed). Convergence is again quite good for 
1000 trajectories. 

Assuming that individual trajectories may have some physical significance, as has been 



suggested 



32| |. it is worth examining a few to see whether their dynamics makes intuitive 



sense. For this moderate barrier case we anticipate asymptotic states which are mixtures of 
left and right. In fact as we see in Fig. 16, where we plot densities for individual trajectories 
at time t = 12, some trajectories do lead to superpositions. However, some also remain 
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strongly localized on the left of the barrier, while others have made a jump to the right. 
Thus, superselection appears to begin to play a role even for moderate barriers. 

Now consider the case of a higher barrier, where e = .1, which is still very modest 
compared to that expected for large chiral molecules. In Fig. 17 we again plot {x\il)t) {'4't\x) 
vs X for individual trajectories, this time at t = 14. For this higher barrier the individual 

FIG. 17: {x\tpt) {ipt\x) vs X for individual trajectories at t = 14, for Example IVA with e = .1 
0.8 I 1 1 1 1 1 1 




-4 -2 2 4 

X 



densities are more strongly chiral, with most strongly localized on one side of the barrier 
or the other. The only state even close to an equal superposition of left and right is the 
dot-dashed curve. These results support the notion of the emergence of a superselection rule 
favoring states localized on one side of the barrier or the other over superpositions. 

In Fig. 18 we follow a single trajectory as it tunnels from one side of the barrier to the 
other. The densities indicated are for times t = (solid), 2 (dashed), 4 (dotted), 6 (dot- 
dashed) and 14 (double-dashed). The tunneling process appears to start almost immediately 
and is effectively over by t = 6. In other cases tunneling had already occurred by t = 2. 
This apparent discrepancy invites further scrutiny in light of the many studies predicting 
well defined tunneling times 
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FIG. 18: {x\^pt) {''Ptlx) vs x for an individual trajectory at various times, for Example IVA with 
e = .1 
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B. Example IVB 

Here we consider a two-level atom immersed in the radiation field of a dielectric band gap. 
This model has recently been studied using approximations to the NMQSD equal ions j2l|. 
In an interaction picture, rotating with the subsystem Hamiltonian, H = and L = ^i'^x — 
i(Jy). The memory function is given by 

a(t,s) = e^(-^)(*-)[Jo(|(t-s))]3 (35) 

where hu is the excitation energy of the atom, a band of allowed frequencies lies between 
A ~ B and A + the gap lies between and A — B, and Jq is the Bessel function of the 
first kind of order 0. We scale time in units oih/\ and energy in units of A. In these units 
we set C(j = 10, 5 = 5 and we will consider two values of A corresponding to uo in the band 
{A = uj) and in the gap {A = uj + B + 1). 
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We fitted to the form (fT^ using nonlinear least squares, where ujj come in pairs 

uj — A — for each and 7^. The resulting parameters are 



u! — A + uj'^ and u;, 



given in Table 1. In Fig. 19 we plot ()35p (solid curve) against (fT3j) (dashed) , which shows 
that the fit is quite satisfactory. 

We chose an initial state \iPq) = |1) which corresponds to an excited 2-level atom (note 
that our convention for labeling the states is the reverse of that in Ref. 3])- We calculated 
(l|pi|l) and {2\pt\2) using 1000 trajectories for u in the band {u = A, sohd and dotted 
curves) and u in the gap {u = A + B + 1, dashed and dot-dashed curves). These quantities 
are plotted in Fig. 20. When u is in the band, emission of a photon occurs, and the two- 
level atom relaxes to its ground state. When u is in the gap the system evolves toward a 
statistical superposition of the two states (i.e. {l\pt\2) = 0) which is weighted toward the 
excited state. It should however be remembered that each individual trajectory represents 
a pure state and so has a density matrix with nonzero off-diagonal elements. In Fig. 21 
we show the density matrix elements (l|pi|l) (solid) and (2|pi|2) (dashed) and the real and 
imaginary parts of (l|pt|2) (dotted and dot-dashed) for an individual trajectory with u in 
the gap. Individual trajectories thus evolve toward a coherent superposition of the excited 
and ground states, with the excited state weighted more heavily. 



27 



FIG. 19: a{t, 0) for Example IVB with lo = A 
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FIG. 21: {l\pt\l), {2\pt\2) and (l|/3t|2) for a single trajectory of Example IVB 




C. Example IVC 



Finally, we consider a three-level system (|1) with energy E'l = in eV, |2) with E2 = 
4.4, and |3) with = 0) representing three electronic states of a Mg''" ion in an ion 
trap, which is driven by a laser resonant with the transition between levels 1 and 2. This 

due to the interesting phenomenon 
of intermittent fluorescence exhibited by the ion. The ion which normally cycles between 
states 1 and 2 occasionally jumps into the dark state 3 giving rise to periods where no 
fluorescence is observed. These stochastic jumps have been observed experimentally^] and 
theoretically^, |^ over a time scale with millisecond resolution. It has been suggested that 
such stochastic jumps might occur on all time scales 5]. Here we show that NMQSD predicts 
jump phenomena on a picosecond timescale consistent with this scenario. 

The effective Hamiltonian in a frame rotating with the Hamiltonian of the isolated ion 
has the form 

iJ=(fi/2)(|l)(2| + |2)(l|) (36) 



where fl is the Rabi frequency of the 1-2 transition. The coupling operators L, for i 
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have the forms 



Li = Ai2|l)(2| 
L2 = Ai3|l)(3| 
L3 = A3i|3)(l| 

L4 = A(|l)(l|-|3)(3|). (37) 



The parameters were chosen as A12 = \Jl /t, A13 = ^R^/r, A31 = ^R^/r where 7 is the 
spontaneous decay rate for the 2-1 transition and -R_ and R+ are the rates out of and into the 

n 

dark state 3, respectively H. It has been shown that -R_ = 8i7^7/9a^ and i?+ = i7^7/18a^ 
where a denotes a Zeeman sphttingj^. Finally, the A operator arises due to interaction 
with the photodetectors. To reproduce the results of Ref. in the Markovian limit we 
set T = /o°° Re a{t)dt . The memory function, which is common to the four noises, was 
calculated with a Debye distribution of frequencies and for a temperature of about .1 K. We 
chose time units 10^^7~^ - roughly a picosecond since 7 = 27r43 MHz[3] - in terms of which 
we set a = 12.1, = 2, A12 = lO'^y^, A31 = ^l/lSn/ (a^^), A13 = ^/8/9n/{a^) and 
A = .22 /sqrtr. We calculated that r = 509 from the 5 term non-linear least squares fit to 
the memory function, although the principle decay occurs over the first 50 time units. The 
dynamics is thus quite strongly non-Markovian. 

In Fig. 22 we show the occupation probability of the 1-2 manifold Kll'i/'i)!^ + |(2|'?/'t)p 
plotted against time. At various times and for various lengths of time this probability 
approaches zero indicating that the ion has jumped to the dark state 3. A more detailed 
analysis of the statistical properties of the dark periods will be presented elsewhere. 



V. DISCUSSION 



Non-Markovian quantum state diffusion (NMQSD), being an exact unraveling of the mas- 
ter equation for an arbitrary subsystem interacting linearly with a boson bath, potentially 
has a wide range of applications in quantum optics. This potential has not been realized to 
any significant extent due to the difficulty of solving - even numerically - the variational- 
differential evolution equation. In fact the development of NMQSD has not led to the 
exact solution of any new problems. In this manuscript we have shown how the variational- 
differential equations (VDEs) of NMQSD can be exactly rewritten as integrodifferential 
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FIG. 22: |(l|V't)P + \{'^\'^h)? vs t for a single trajectory of Example IVC 
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equations which can in turn be rewritten as ordinary or partial differential equations. We 
illustrated application of the new NMQSD equations by solving a number of problems for 
which exact solutions were already known. Both linear and nonlinear versions of NMQSD 
were studied. We found that both versions worked well and yielded high accuracy solutions. 
Finally, we applied the method to three previously unsolvable problems (tunneling in a dou- 
ble well, two-level atom in a photonic band gap, and intermittent fluorescence in a driven 
three-level ion) to show that interesting work can be done with the reformulated theory. We 
anticipate that other interesting applications of NMQSD will emerge now that the equations 
can be solved in a systematic fashion. 

The authors acknowledge the support of the Natural Sciences and Engineering Research 
Council of Canada. 
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